#'
#'
#' An interface of point pattern reconstruction software writen by Prof. Thorsten Wiegand.
#' 
#' @usage
#' ppreconstruction(x,smst=list(),nrc=99,intmap=TRUE,
#'  keepfiles=TRUE,dimyx=c(10,10),sigma=5,read_rec=TRUE)
#' 
#' @param x a population or community object
#' @param smst a list of settings about summary statistics used in reconstruction. in current version, it only use the default summary staistics in the softeware.
#' @param nrc number of populations to be reconstructed
#' @param intmap the map of point intensity in reconstruct a heterogeneous point pattern
#' @param keepfiles whether to keep middle files created by Thorsten's software
#' @param dimyx the grid size of the individual intensity map
#' @param sigma the smooth width in the kernel smooth estimation of the intensity map.
#' @param read_rec a logical flag to read the reconstruction file. you can read the reconstructed files after the reconstruction process for effictive use of physical memnry.
#' 
#' @return
#' a reconstructed population
#' 
#' @note
#' Please note that this only works on the windows plantform. please pay atterntion for the possible warnings
#' 
#' @examples
#' data(BCI)
#' pop1=as.population(subset(BCI,BCI$species=="DES2PA"))
#' attr(pop1,"species")="as*dg**s-U--s" #testing for illegan file name
#' rec_pops=ppreconstruction(pop1,nrc=1,dimyx=c(100,100),sigma=0.5,keepfiles=FALSE)
#' plot(pop1)
#' plot(rec_pops[[1]])
#' 
#' com=remove_rare_species(BCI,1600)
#' rec_coms=ppreconstruction(com,nrc=1)
#' plot(com)
#' plot(rec_coms[[1]])
#' 
#' 
#'spab=species_abundance(BCI)
#'testsp=names(spab)[spab>=56 & spab<57]
#'subbci=get_populations(BCI,testsp)
#'re_com=ppreconstruction(subbci,nrc=2)
#'
#'
#'#test the type I error of this method
#'#step1: generate a community with ten independent species. Each species is nonrandomly distributed 
#'# The best way to do it is torrus translation
#'data(BCI)
#'BCI=remove_rare_species(BCI)
#'spab=unique(sort(species_abundance(BCI)))
#'spchosed=names(spab)[c(seq(1,length(spab),12),length(spab))]
#'bci=get_populations(BCI,spchosed)
#'
#'#pvalues=list()
#'
#'for(i in 1:199){
#'
#'bci_obs=rtorus(bci)
#'
#'#step2: call  test_pairwise_association_community to get a pvalue
#'#re=test_pairwise_association_community(bci_obs,nsim=199,r=0:25,dimyx=c(20,20),sigma=10)
#'#pvalues[[i]]=lapply(re,function(x) attr(x,"pvalues"))
#'re_com=ppreconstruction(bci_obs,nrc=199,keepfiles=FALSE,dimyx=c(20,20),sigma=10,read_rec=TRUE)
#'save(re_com,file=paste("simu",i,".RData"))
#'}
#'
#'#call test_pariwise_association to calcualte type I error.
#'
#'#step3: repeat step1 and step2 for 199 times, then we can calcualte the type I error of this point pattern reconstruction method
#'
#'all_rec=list()
#'
#'
#'
#'
#'

ppreconstruction=function(x, ...){
  UseMethod("ppreconstruction")
}

ppreconstruction.community=function(com,smst=list(n_max=40000,min_error=0.005,ar_ratio=0.01),
                                    nrc=99,intmap=TRUE,keepfiles=TRUE,dimyx=c(10,10),sigma=5,read_rec=TRUE){
  pops=split(com)
  S=total_richness(com)
  rec_pops=list()
  cat(paste("Total",S,"species,","finished species:"))
  for(i in 1:S){
    rec_pops[[i]]=ppreconstruction(pops[[i]],smst,nrc,intmap,keepfiles,dimyx,sigma)
    cat(paste(i,",",sep=""))
  }
  rec_coms=list()
  if(read_rec){
    for(i in 1:nrc){
      rec_coms[[i]]=lmerge(lapply(rec_pops,function(x) x[[i]]))
    }
  }
  
  return(rec_coms)
}

ppreconstruction.population=function(pop,smst=list(n_max=40000,min_error=0.005,ar_ratio=0.01),
                                     nrc=99,intmap=TRUE,keepfiles=TRUE,dimyx=c(10,10),sigma=5,read_rec=TRUE){
  if(.Platform$OS.type!="windows")
    stop("Can't run the pattern reconstruction software under non-windows OS")
  
  spname_orig=attr(pop,"species")
  #remove illegal characters in the spname
  spname=gsub("([*])","",spname_orig)
  spname=gsub("([-])","",spname)
  
  #create a work dir
  tempdir=paste("./pprTemp",spname,sep="_")
  isdone=dir.create(tempdir)
  setwd(tempdir)
  
  pre_filename="a"
  
  #find the already done jobs and start job id
  startid=findstartid(pre_filename,nrc)
  nleft=nrc-startid+1
  
  if(nleft>0){
    
    #get the package library path
    pkgpath=find_package_path("scp")
    
    #step 1, prepare a setting file for reconstruction
    file.copy(paste(pkgpath,"/exec/default.rec",sep=""),"./default.rec")
    modify_setting_file(pop,smst)
    
    #copy the exe program into current wd
    file.copy(paste(pkgpath,"/exec/PatternReconstructionMay2015.exe",sep=""),
              "./PatternReconstructionMay2015.exe")
    
    #step 2, prepare a data file
    spintensity=density(pop_to_ppp(pop),dimyx=dimyx,sigma=sigma)
    write_prerec_file(pop,covr=spintensity,filename=pre_filename)
    
    #step 3, start the program in batch model
    #it seem there are some format error in the prepared file, the system command is no problem
    try(system(paste("PatternReconstructionMay2015.exe",
                     "default.rec",
                     paste(pre_filename,".dat",sep=""), 
                     paste("nr*",startid," rep*",nleft,sep=""),
                     paste("int_",pre_filename,".int",sep=""))))
  }
  
  #step 4, read the reconstructed patterns
  if(read_rec){
    pops=try(read_multi_constructions(paste("int_",pre_filename,sep=""),spname,plotdim(pop),nrc))
    pops=lapply(pops,function(x) {attr(x,"species")=spname_orig;return(x)})
  }else{
    pops=c()
  }
  
  
  #step 5, remove middle files if required
  if(!keepfiles){
    file.remove(dir())
    setwd("../")
    shell(paste("rd",substr(tempdir,3,nchar(tempdir))))    
  }else{
    file.remove(c("PatternReconstructionMay2015.exe","default.rec","a.dat","int_a.int"))
    setwd("../")
  }
  
  #step 6, return the results
  return(pops)
}

write_config_file=function(smst,pre_filename){
  
}

#it will find the maximum done id and return a next integer
findstartid=function(pre_filename,nsim){
  prec_rec=paste("rec_int_",pre_filename,"_",sep="")
  for(i in nsim:0){
    filename=paste(prec_rec,i,".dat",sep="")
    isdone=file.exists(filename)
    if(isdone){
      return(i+1)
    }
  }
  return(i+1)
}

# pop a populatin object
# covr (optional) a intensity map, e.g. intensity of all points or effect map of covraiates
# filename the name of the point coordination file. it should be a "*.dat" file. if covr is provided, the file of covr will be "cov_*.dat". path of the intensity file will be as same as the give dat file.
# LINUX whether it is on a linux system
write_prerec_file=function(pop,covr=NULL,filename,LINUX=FALSE,...){
  #setup the line break for windows user
  if(LINUX)
    eol="\r\n"
  else
    eol="\n"
  plotdim=plotdim(pop)
  plotinf=matrix(c(0,plotdim[1],0,plotdim[2],total_abundance(pop)),nrow=1,ncol=5)
  write.table(plotinf,paste(filename,".dat",sep=""),row.names=FALSE,col.names=FALSE,eol=eol,...)
  points=data.frame(pop$x,pop$y)
  write.table(points,paste(filename,".dat",sep=""),row.names=FALSE,col.names=FALSE,eol=eol,
              append=TRUE,quote=FALSE,...)
  if(!is.null(covr)){
    temp=strsplit(filename,split="/")[[1]]
    temp=temp[length(temp)]
    pre=strsplit(filename,split=temp)[[1]]
    temp=strsplit(temp,split=".dat")
    filename.covr=paste(pre,"int_",temp,".int",sep="")
    ngrid=covr$dim
    grid_size=covr$xstep
    plotinf2=c(1,ngrid[2],1,ngrid[1],prod(ngrid),grid_size)
    plotinf2=matrix(plotinf2,nrow=1)
    xy=expand.grid(y=1:ngrid[1],x=1:ngrid[2])
    xy=xy[,c(2,1)]
    #normalize the covraite map
    covr$v=covr$v/max(covr$v,na.rm=T)
    density=cbind(xy,rep(1,dim(xy)[1]),round(as.numeric(covr$v),6))
    write.table(plotinf2,filename.covr,row.names=FALSE,col.names=FALSE,eol=eol,...)
    write.table(density,filename.covr,row.names=FALSE,col.names=FALSE,eol=eol,
                append=TRUE,quote=FALSE,...)
  }
}

# filename the path and filename of the reconstructed file
# return value a conrespondence ppp object
read_rec_file=function(filename){
  title=scan(filename,nlines=1,sep=" ",quiet=TRUE)
  title=title[!is.na(title)]
  window=owin(c(title[1],title[2]),c(title[3],title[4]))
  nlines=title[length(title)]
  xy=scan(filename,nlines=1+nlines,sep=" ",quiet=TRUE)
  xy=xy[!is.na(xy)]
  xy=xy[-c(1:length(title))]
  dim(xy)=c(4,length(xy)/4)
  xy=t(xy)[,1:2]
  data.ppp=ppp(x=xy[,1],y=xy[,2],window=window)
  return(data.ppp)
}


# prefilename the prestring of all reconstructed files becomes to a simple plot
# nsim number of repeat for each species
read_multi_constructions=function(prefilename="rec_bci_rec",species,plotdim,nsim=99){
  
  null_communities=list()
  for(j in 1:nsim){
    rec_filename=paste("rec_",prefilename,"_",j,".dat",sep="")
    tp=read_rec_file(rec_filename)
    null_community=population(species,x=tp$x,y=tp$y,plotdim,valid_check = FALSE)
    null_communities[[j]]=null_community
  }
  return(null_communities)  
}


#help to find the location of the installed scp package
find_package_path=function(pname){
  #try the default search place, R_HOME
  ppath=find.package(pname,quiet = TRUE)
  if(length(ppath)==0){
    longversion=getRversion()
    newsearchpath=paste("~/R/win-library/",longversion,sep="")
    ppath=find.package(pname,newsearchpath,quiet = TRUE)
    if(length(ppath)==0){
      shortversion=substr(longversion,1,3)
      newsearchpath=paste("~/R/win-library/",shortversion,"/",sep="")
      ppath=find.package(pname,newsearchpath,quiet = TRUE)
    }
    if(length(ppath)==0){
      stop("can not find the installed scp packages")
    }
  }
  return(ppath)
}

modify_setting_file=function(pop,smst=list(n_max=40000,min_error=0.005,ar_ratio=0.01)){
  con=file("default.rec")
  sets=readLines(con)
  sets[6]=paste("You calculated summary statistics up to ",round(min(plotdim(pop))*0.4,0)," data units",sep="")
  writeLines(sets,con)
}
# 
# write_setting_file=function(pop,smst=list(n_max=40000,min_error=0.005,ar_ratio=0.01)){
#   words="\nThis file gives you the settings and the results of the pattern reconstruction of file a.dat\n"
#   words=paste(words,"You stopped reconstruction after ",smst$n_max," iterations\n",sep="")
#   words=paste(words,"You stopped if the total energy dropped below ",smst$min_error,"\n",sep="")
#   words=paste(words,"You stopped if the ratio of accepted/proposed points dropped below ",smst$ar_ratio,"\n",sep="")
#   words=paste(words,"You calculated summary statistics up to ",round(min(plotdim(pop))*0.4,0)," data units\n",sep="")
#   
# words=paste(words,"You used bins of 1 data units
# You used the option for strong aggregation where in 15  of all cases a new tentative point is placed within radius 3  of an already placed point
# You did not use a torus shift as start pattern
# You used every simulation a different sequence of random numbers
# you conducted homogeneous pattern reconstruction without intensity function
# You used 45 test points in one direction for calculating Hs(r)
# You used a grid of test points for estimation of Hs(r)
# You used no edge correction for the nn statistics
# You used estimator method 2 for the pair correlation function
# 
# You used the following test statistics:
# Pkr(k,r=5)
# Pkr(k,r=15)
# Pkr(k,r=45)
# r*g(r)
# the NN distance distribution D(r)
# the K-function K(r)
# the spherical contact distribution Hs(r)
# the distribution D1(r) to 1st neighbor
# the distribution D2(r) to 2st neighbor
# the distribution D4(r) to 4st neighbor
# the distribution D6(r) to 6st neighbor
# the distribution D8(r) to 8st neighbor
# the distribution D12(r) to 12st neighbor
# the distribution D16(r) to 16st neighbor
# the distribution D20(r) to 20st neighbor
# the distribution D25(r) to 25st neighbor
# the distribution D50(r) to 50st neighbor",sep="")
#   
#   write.table(words,file="default.rec",quote=FALSE,row.names = FALSE,col.names = FALSE)
#   return()
# }